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The dissipative collision of two identical viscoelastic disks is studied. By using a known law for 
the elastic part of the interaction force and the viscoelastic damping model an analytical solution 
for the coefficient of restitution shall be given. The coefficient of restitution depends significantly 
on the impact velocity. It approaches one for small velocities and decreases for increasing velocities. 



I. INTRODUCTION 

The coefficient of restitution is an important means to characterize the damping properties of granular particles. It 
is defined as 

with g being the absolute value of the normal component of the relative velocity before and g' the value after the 
collision. In most analytical and numerical studies of the behaviour of dilute or fluidized granular systems a constant 
coefficient of restitution was used (see, e.g., 0, S, 0, II 0, 1 and many more). In experiments on (3d) spheres, however, 
it has been found that this coefficient is not a constant but depends significantly on impact velocity 0- For spheres 
there is a theory based on first principles which predicts this velocity-dependence @, Q . The aim of the present paper 
is to use the methods successfully applied to spheres for the description of the collisional properties of disks. 

If we assume particles whose conservative interaction part is linear with respect to the mutual compression and 
whose dissipative interaction part is linear with respect to the relative velocity the coefficient of restitution would 
indeed be constant (Io| . This can be very useful if one wishes to compare results from Molecular Dynamics (or Discrete 
Elements) simulations with the predictions from theory based on e =const. Apart from these practical considerations 
there are authors who use disks as a real world examples of particles interacting via a linear force law (see, e.g., 
[HI 03 )• However, it is known for several decades that the contact force law for colliding disks differs significantly 
from a linear law [l3j: 
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Here £ = 2R — \r\ — T2\ is the compression of the particles of radius R at positions ri/2- The material properties are 
characterized by the Young modulus Y and the Poisson ratio v. One should note that in the present context "disk" 
means infinetely long cylinder. Above expression is, thus, only approximately valid for disks of finite thickness (or 
cylinders of finite length). For too thin disks the (implicit) force law ([2|) may become wrong. Eq. ([2]), can be solved 
for the contact force F e i yielding 



The function Wq is the zeroth Lambert-IT-function. It has two real and infinitely many complex branches. One of 
them vanishes linearly as the argument vanishes, the others diverge as the logarithm of the modulus of the argument. 
The first (vanishing) branch is irrelevant since it would yield a finite force for infinitesimal deformation. The relevant 
details of Wo are summarized in appendix [A] The force law (J3j) was used to study the coefficient of restitution of 
colliding disks where energy was dissipated via excitation of internal degrees of freedom of the disks, i.e., by excitation 
of vibrational modes [l4[ • This dissipation mechanism will be neglected here. 
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In this article we assume viscoelastic damping of the deformed particle material. For this damping mechanism 
the material of contacting particles is assumed to be locally linear, i.e. the stress tensor depends linearly on both 
the deformation tensor and the deformation rate tensor. If the impact velocity is small enough this assumption is 
reasonable. For higher velocities other dissipation mechanisms, like plastic deformation or the already mentioned 
excitation of vibrational modes have to be taken into account. The authors of [1 51 ] used the viscoelastic model to 
study the contact of spheres and found that the dissipative force component Fdi ss is related to the conservative force 
F e \ via 

F diss = Ai^, (4) 

where (in 3d) F e \ is the Hertz contact force. Subsequently the above expression Eq. (U) was derived using a crystal 
mechanical approach [l6j . The parameter A is proportional to the viscous relaxation time of the particles and depends 
only on the material constants and the geometry, such as the radius of the particles. Due to the general nature of the 
argumentation from [l5T ] it can be applied to the colliding disks problem too. Inserting the force law, Eq. ([3]), into 
Eq. g]) yields 

Fd\ ss — -, „ +1 x • (5) 

1 + W 

With Eq. ([4]) the equation of motion reads 



4i? 



m cS t + F cl + A£—f = (6) 

m = o (?) 

e(0) - g (8) 

where repulsive forces are positive. We introduce rescaled variables for length, force and time 

* = ST « (9 » 

F e i = F e i (10) 

t 11 

m 

and define the scaled velocity v and the scaled damping parameter a by 

v = — = g J (12) 

dr y 4R V K ' 

a = aM- . (13) 
V m 

Introducing Eqs. ©-(13]) into ©-© we obtain 

x + F cl + ax^- = (14) 
dx 

x(0) = (15) 

i(0) = v , (16) 

There is a simple interpretation of v and a. We can rewrite m — p ttR 2 (with p being the material density) and find 
v ~ 9y p/Y . In a rough approximation one can estimate the speed of sound by c 

9 

v ~ — . 

c 



(17) 



By a similar argumentation we find that a ~ Ac/R, i.e., a can be estimated by the viscous relaxation time (which is 
the meaning of the damping parameter A) divided by the time a sound wave needs to travel through the cross-section 
of the particle. Obviously both v and a must be small compared to unity. 
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The implicite elastic force law in the rescaled variables does not contain any dependence on material parameters 

x = -F el lnF sl . (18) 

The explicit solution of Eq. (fT8|) is 



F, 



el 



Wb (-as) 



(19) 



and can alternatively be obtained substituting the normal variables by the rescaled ones in Eq. ((31). The rescaled 
force as a function of the rescaled compression is given in Fig. [1] 
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FIG. 1: The rescaled elastic force F e \(x) as a function of the rescaled compression x. 

The dissipative part of the interaction force (see Eq. ([4}) is 

dF„] x x 



dx 



1 + lnR, 



l + Wo(-x) 



(20) 



In a straightforward approach the coefficient of restitution can be found by integrating Newton's equation of motion 
and determining the trajectory of the particles x(t), the duration of the collision r c and, eventually, e = —x(t c )/v. 
Here, however, we use a different approach due to mathematical difficulties resulting from the complicated form of 
the force law. 

In the next section [TT] an approach which doesn't rely on the explicit knowledge of the trajectory is presented. The 
method allows to compute the dominant part of the energy loss during the collision. As a further advantage we do 
not need to know explicit expressions for the damping force, we only need to know its principal structure ([3]), which 
is especially useful for more complicated force laws. In section IIIII a method to further improve the obtained result 
will be presented. 



II. ENERGY BALANCE 



The coefficient of restitution and the total energy loss during a collision are related by 

, AE 
1 - e 2 = 2^- . 



(21) 



The value of AE up to linear order in the damping parameter a can be found without explicitly knowing the trajectory 
of the particle. This procedure was successfully employed before to compute the coefficient of restitution for colliding 
spheres [H [t| [13, EH • We introduce the interaction potential <!> and rewrite the equation of motion (Q2 



_d_ 
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■ dF i 
dx 



(22) 



4 



introducing the potential <£> via 

= die ^ ' (23) 

Note that the custumary negative sign is missing in the above definition as we count repulsive forces positive. The 
potential energy shall vanish for vanishing deformation. The term in brackets on the left hand side of Eq. (f2"2")l is the 
total energy, 

E = Y + <P (24) 

(with initial value Eq = v 2 /2). Thus we can rewrite Eq. (|22[) and obtain the energy balance equation 

^ = ~ax 2< ^^- = ~F diBS i . (25) 
dr ax 

For the total energy loss during the collision we thus find 

AE = -a[dT± 2 ^, (26) 
J dx 
o 

where t c is the duration of the collision. To actually evaluate this integral we would now need expressions for the 
trajectory of the particles. However, for small damping the trajectory will be close to the undamped trajectory xq, 
i.e. x(t) = xq(t) + axi(r) + 0(a 2 ), where the first order correction term x\ is well behaved for all practical collision 
problems. Inserting this ansatz into Eq. ([26]) we find that we can estimate the energy loss in linear order of a by 
evaluating the integral along the undamped trajectory instead of the correct damped one. We thus find 

AE = -a J dr x 2 ^- + O (a 2 ) , (27) 
o 

where t,? is the duration of the undamped collision. The force F e \ has to be evaluated at the positions of the undamped 
trajectory. This approximation allows us to determine AE without explicitcly knowing the trajectory. We split the 
integral into two parts 

AE = -a [ drx 2 ^-a [ drx 2 ^ + O (a 2 ) , (28) 
J dx J da; 

r°/2 

where r^/2 is the time of maximum compression, where the particles are momentarily at rest and start to seperate 
again. Since the trajectory is symmetrical with respect to this time the second part of the integral is exactly equal to 
the first. 

AE = -2a [ drxl^ + Oia 2 ) (29) 
J dx 
o 

Since in the interval (0,r°/2) the force is a monotonous function of time we can replace the time integral by an 
integral over the (elastic) contact force, i.e., 

drx V 1 = dFel ■ (30) 
dx 

Furthermore we express the velocity of the undamped problem in terms of the potential energy, which gives 



±o = y / 2(E-$) = Vv 2 - 2$ , 



(31) 
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with v being the (scaled) impact velocity. Therefore, 



-fmax 

AE=-2a J (LFdvV -2$ + C(a 2 ) 



(32) 



We express the potential energy in terms of the force, i.e., $ = $(</>)■ 

- d$ d$ dF i 

^el — -3— — — ; — 

da; dFci da; 



dF el 

From Eq. (JTSJl we find dF ol /dx = -1/ (l + lnF e i) 

d$ 



and thus 



dF nl 



= -F el l+lni^ e i 



With $(i = 0)=0we find 



Therefore 



$ = -^f (l + 21nF e i 



AE = -2av 



/ M\ 1 



F e ^l + 21nFei 
2^2 



+ 0(a 2 ) 



(33) 
(34) 



(35) 



(36) 



(37) 



The maximum elastic force -F max is related to the impact velocity by inserting the energy conservation v 2 = 2$ into 
Eq. 436} 



2v 2 = ~F^Jl + 2\nF, 



(l + 21nF max ) 



(38) 



We will now use this expression to replace v in the integral. The actual solution F max (v) will be given at a later stage 
of the computation. 



AE = -2av 



F 2 (l + 2lnF cl 
^ ^ ax (l + 21n^ max ) 



+ 0{a 2 ) 



Substituting F c \ = -Fmax y with y from the interval (0, 1) we find 



AE = -2av4> mgx J dy 
\ 



1 - 



i/ 2 (l + 21nF nuut + 21ny 



1 + 21x1^ 



+ 0(a 2 ) 



-2av(f) max / dy 



\ 



l-y' 



2y 2 In y 



1 + 21nKr 



+ 0{a 2 ) 



1 

= -2av4> max J dyy/l- y 2 ^ 



1 - 



2y 2 In y 



1 + 2111^ 



0(a 2 



(39) 



(40) 

(41) 

(42) 
(43) 
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In the interval y £ (0, 1) the relation 



2y 2 In y 



< 1 . 



(44) 



holds. For sufficiently small impact velocity v we can, therefore, expand the second term in the integrand into a power 
series with respect to the small term 2y 2 lny/ [(l — y 2 ) (1 + 2m<^ max )l . With a k being the Taylor coefficients of the 



expansion of y/1 + x around i = 0we find 

oo \ 

AE = -2awF max ^(-l) fc afc / dy y/l - y 2 

I n J 



2y 2 In y 



= -2avF max ^ 



{-l) k a k 



(1-y 2 ) (l + 21nF„ 
k 



+ 0{a 2 ) 



-2avF max ^2 



k=0 ( 1 + 21nF max ; n 
(-l) k c k 



k=0 ^ 



1 + 21nF n: 



dyy/i - V 2 

u+0{a 2 ) , 



V hit/ 
i-y 2 



+ 0{a 2 ) 



(45) 
(46) 
(47) 



with 



Ck = a k 



dyVl 



2y 2 \ny 



i-y 2 



(48) 



For ci and c-i there exist analytical expressions, the higher coefficients have to be computed numerically. The first 
values of c k are given in Table U at the end of next section. We now have to find an expression for F ma x in terms of 
v. The defining equation ([55]) can be solved in closed form using again Lambert's IF-function: 



2v 2 



l + 21nF max = W {-2ev 2 ) 



(49) 
(50) 



with the Euler constant e. Note that the function Wq is negative for negative arguments. The final expression for the 
energy loss reads 



AE = -2av 2 . 



-W (~2ev 2 ) ' [W (-2ev 2 )] k 



k=0 
C k 



ifc+1 



+ 0(a z ) . 



kTo [~W (~2ev 2 )Y 

The coefficient of restitution can now be determined with the same accuracy, i.e., linearly in a. We have 



e = \ 1 



AE 
~E~ 



2AE 



1 + 



AE 



(51) 
(52) 

(53) 
(54) 
(55) 



Inserting Eq. 
parameter a. 



we arrive at the final expression for the restitution coefficient up to linear order in the damping 



1 -V8aJ2 



^ [-W (-2ev*)] 



(56) 
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For very small velocities this expression can be approximated as 



1-e 1 (57) 

^/-W (-2ev*) 

~ ; 1 (58) 
v /ln(l/2eu 2 ) 

in accordance to the previously published result of Hayakawa and Kuninaka [l8j . However, for the approximation 
(|58|) to be useful the expression ln(l/2eu 2 ) has to approximate the Lambert W^-function to a reasonable accuracy 5. 
To this end we have to require 

1 >\(*J\' (59, 



y/\n(l/2ev 2 ) S yy/]n(l/2ev 2 ) 
1 
6 

v < e 1 / 25 " ' 5 (61) 



ln(l/2eu 2 ) > j (60) 



The critical velocity becomes exponentially small for increasing accuracy requirements. This severely limits its appli- 
cability. The full solution (to linear order in a) Eq. (|56|) does not suffer from this limitation. 

The energy balance method unfortunately doesn't allow the calculation of higher order contributions to the energy 
loss, since this would require detailed knowledge of the trajectory. However, there is a method to calculate the 
contribution proportional to a 2 which will be described in the next section. 



III. CONSISTENCY METHOD 



In the previous section we have seen that the dominant term of the restitution coefficient is linear in the damping 
coefficient a. With the help of a novel method it shall be attempted to improve the result to the next order in a. 
It will be shown that for any viscoelastic problem of form the second order contribution to the coefficient of 
restitution is completely determined by the first order contribution. We assume that the next largest contribution is 
quadratic in a, or generally that the restitution coefficient e can be expressed as a power series in a: 



e(v) = J2<* k fk(v) (62) 



k=0 

Note that the functions do not depend on a. The zeroth function /q is 1 since the restitution coefficient has to be 
one for undamped collisions (a = 0). After the collision the relative velocity of the two particles will be 

v' = e(v)v . (63) 

Now we perform a time inversion, i.e. we start with the velocity v' at time t c and let the time run backwards until 
at time t = we arrive again at the initial velocity of the direct collision v. The equation of motion of this inverse 
problem reads 

x + 4>~ a-f- = 64) 

ax 

Comparing with the equation of motion of the direct (forward time) collision (|14p one notes that the only difference is 
the sign of the damping coefficient a. In the inverse problem we have a negative damping —a so, in accordance with 
intuition, the inverse collision is an accelerated one. We now know the restitution coefficient of the inverse problem 



r (v) = J2(-a) k Mv) > (65) 



k=0 



again the only difference being the opposite sign of a. As said before, starting the inverse collision at velocity v' we 
must arrive at the velocity v: 

v = e inv (v') v' . (66) 
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Inserting v' from Eq. (|63f into t|66|) we obtain the consistency equation 

£inv [ve(v)] e(v) = 1 . 



(67) 



Eq. (|67|) has to be fulfilled for all restitution laws of form Eq. ([62]) regardless of the internal contact mechanism. 
Of course, different contact laws yield different velocity dependences for e(v) and e±ny(v), however, the consistency 
requirement Eq. (|6T[) in it's general form remains unchanged. Eq. (|67p allows us to calculate the second order 
correction to the first order result Eq. (|56|) . We insert Eqs. (|62j) and (|66|) into Eq. ([67]) using /q = 1 and 



/ fe [ve(v)] = f k 



v + ^a k vf k (v) 



k=l 



/*(«)+ -^rWE a ^w + 2d7 



fc=i 



\fc=i 



and find 



1 = 1 - of 



(vh{v)^{v) + f?(v)^j + 2a 2 f 2 (v) + 0(a 3 ) 



(68) 
(69) 

(70) 



All terms on the right hand side of linear or higher order in a have to be zero which gives us the final expression for 
the second order term / 2 



vfifl + fi 



(71) 



where f 1 means derivative with respect to the velocity v. When computing higher orders one notes that, unfortunately, 
terms of odd order in a vanish identically. We cannot, therefore, compute terms of e of order higher than a 2 since 
half of the necessary equations are missing. From Eq. (|56p we have 



/i(«) = -Vsg. 



^ {-W (-2ev2)] k+1/2 ■ 
In order to calculate the first derivative of this expression with respect to v we use 

dW (x) _ Wo Or) 



and find 



dx x(l + W {x)) 
32 ^ c k (k + |) 



E 



v [1 + W (-2ev 2 )} ^ [- Wo (-2ev 2 )} k+1/2 
With this expression we can calculate vf\f[ 
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(72) 



(73) 



(74) 



vfifl 



1 + W (-2ev 2 



T tzrY(k-i + -) 

to l-W (-2ev*)] k+1 U V V 



CiCk- 



(75) 



Finally f 2 reads 



A 2 = 8E' 



fe=0 



-W (-2ew 2 )] fc+ 



tE c < 



/2 can be expressed in the form 



E 



1 + W (-2ev2) ^ [-W r (-2e« 2 )] fc ' 



(76) 



(77) 
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k 


Cfe 


d k 




1 

2 
3 


7T 

4 

|(l-ln4) 

-0.02237433 
-0.0076646 


16 

-g(l-ln4) 

0.250419 
0.029518 



TABLE I: First values for c k and d k (explanation given in text). 



where the d& are purely numerical constants wich can be computed from the Cfc. 



d 



-2^(fc - i - l)ciCk-i-i - ^^CiCk-i for k > 



(78) 
(79) 



i=0 



i=0 



The first values of Ck and dj, are given in Table U For k = and fc = 1 analytical expressions for both Ck and dfc can 
be found, the other constants have to be determined numerically. 

We now have the final expression for the coefficient of normal restitution up to second order in a: 



e = 1 - y/8a^T- 



Ck 



4a 2 



-W (-2ev 2 )] k+1/2 1 + W (-2ev 2 ) ^ [- Wo (- 2 ev 2 )} 



E 



(80) 



To check the analytical result the equation of motion Eq. (fl~4"|) has been integrated numerically for different values 
of the scaled impact velocity v. The result for a = 0.1 is shown in Fig. [2] (full line). Both analytical and numerical 
solution agree very well. 
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FIG. 2: Coefficient of restitution of two colliding disks for various values of a. The full lines show the result of the numerical 
integration of Eq. (|14p . The dashed lines are the analytical solution up to second order in a (see Eq. I|80p ). From top to botton 
the values of a are 0.02, 0.05, 0.1 and 0.2. To demonstrate the effect of the second order contribution in a the first order result 
is shown with dash-dotted line. One notes a large difference to the numerical result. 



IV. CONCLUSIONS 



In the present article the velocity dependence of the coefficient of normal restitution of colliding identical disks was 
derived. It turns out that it shows a significant dependence on velocity. It approaches one (no damping) for small 
velocities and decreases for increasing velocities. Consequently, for freely cooling systems it is problematic to use 
colliding disks as an example for particles with constant restitution coefficient. However, for driven systems the 
assumption of constant restitution may remain justified. For the computation of the coefficient of restitution it is 
not necessary to know exact details about the trajectories of the particles during the collision but instead a simple 
energy balance method together with a consistency consideration is sufficient to derive an analytical expression which 
is correct up to the second order of the damping parameter. It can be expected that the kinetic properties of gases of 
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colliding disks differ significantly from the properties of gases of particles colliding with constant restitution coefficient. 
The analytic results for e(v) have been compared with results obtained by numerically integrating the equation of 
motion of the collision problem. Both solutions are in good agreement. 

The author wants to thank Thorsten Poschel and Nikolai Brilliantov for helpful discussions. 

APPENDIX A: LAMBERT ^-FUNCTION 

The Lambert VF-function of argument x is implicitely defined by the following equation [l^ | 

We w = x . (Al) 

The defining equation has one real solution for positive values of x but two for negative x. Fig. [3] shows a plot of this 
function for the relevant case x < 0. 




FIG. 3: Lambert VK-function for negative argument. The dashed line is the (irrelevant) branch which approaches zero for 
x — > 0, the full line is the branch which displays the required behaviour. 



One solution approaches zero like W(x) — > x for x — ► 0, the other approaches — oo. It is clear that the first solution 
(W(x) — > x) is unphysical since it would yield a force — > 1 as X — » 0, i.e., a nonvanishing contact force for zero 
compression. To avoid confusion between the two branches the relevant solution will be called Wq instead of W . For 
small negative values of x the function Wq can be approximated by 

W « -ln(--J -lnln(-i J (A2) 



It shall be mentioned that there is no real Wo for x < — e _1 (see again Fig. [3|. This, however, is of no impact to our 
solution since the we are restricted to small arguments both in the force law Eq. (JTUJ) as well as in the final solution 
Eq. ([50)1 . More specifically, an argument x = — e _1 in Eq. (fH)|) would correspond to compressions comparable to the 
particle radius, in Eq. (|80p it would correspond to velocities comparable to the speed of sound in the material, both 
cases are clearly beyond the limits of viscoelastic approximation. 

The first derivative of the Lambert IF-function with respect to its argument can be calculated by differentiating 
both sides of Eq. (|A"Ij) : 



dW 
dx 



e w (l + W) = 1 
dW 



dx e w (l + W) 



or using the defining equation again 



. w 



W 

dW W 



dx x(l + W) 



(A3) 
(A4) 

(A5) 
(A6) 
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These expressions hold for both the real VK-functions (W and Wo). 
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